# The following script was adopted from Thiele, Jan C., Kurth, Winfried and Grimm, Volker (2014) 
# 'Facilitating Parameter Estimation and Sensitivity Analysis of Agent-Based Models: 
# A Cookbook Using NetLogo and 'R'' Journal of Artificial Societies and Social Simulation 17 (3) 11 
# <http://jasss.soc.surrey.ac.uk/17/3/11.html>.


library(RNetLogo)

# an R random seed (for beeing reproducible)
set.seed(-402223867)

# NetLogo path; NetLogo model path; simulation function path

nl.path     <- "C:/Program Files (x86)/NetLogo 5.2.0"
model.path  <- "C:/Users/.../ABSAM_R.nlogo"
simfun.path <- "C:/Users/.../sim_fun1_all.R"

# min. and max. values for the each of parameters      
  input.values <- list(
    'shock'  = list(min = 0.01, max = 0.04),
    'minimum-wage' = list(min = 0.8 , max = 1.7),
    'efficiency'    = list(min = 0.10,  max = 0.30),
    'beta'          = list(min = 0.4, max = 0.6),
    'growth-rate'    = list(min = 0.002,  max = 0.03),
    'benefits'      = list(min = 0.4, max = 1.2)
                    )

# give details for the sampling design of the Morris method
morris.design <- list(type = "oat", levels = 6, grid.jump = 3)

# give the number of repitions of morris screening (argument r in function morris)
no.repitions <- 30

# names of output values
output.names <- c("umean.criterion","theta.criterion","ltu.criterion", "wages.criterion", "ltu.wages.criterion", 
                  "duration.criterion", "ltu.duration.criterion")

# how many repetitions for each parameter set should be run

no.repeated.sim <- 20

plot.on.screen <- TRUE

# should R report the progress
trace.progress = TRUE

# initialize NetLogo
nl.LMA2.0R <- "nl.LMA2.0R"
NLStart(nl.path, gui=TRUE, nl.obj=nl.LMA2.0R)
NLLoadModel(model.path,nl.obj=nl.LMA2.0R)


# load the code of the simulation function (name: simulate)
source(file=simfun.path)

require(sensitivity)

# variable used for progress tracing
already.processed <- 0  
  	
# get names of parameters
input.names = names(input.values)

# calculate number of iterations
iter.length <- no.repitions * (length(input.values)+1)

# get the min and max values of the input factor ranges
mins <- sapply(seq(1,length(input.values)), function(i) {
		input.values[[i]]$min})
maxs <- sapply(seq(1,length(input.values)), function(i) {
		input.values[[i]]$max})

# create input sets
mo <- morris(model = NULL, factors = input.names, r = no.repitions, design = morris.design,
            binf = mins, bsup = maxs, scale=TRUE)

# simulate for all input sets
# get results of all evalulation criteria as matrix                                                    
sim.results.morris <- apply(mo$X, 1, function(x) {simulate(param.set=x,                            
                              no.repeated.sim=no.repeated.sim, 
                              nl.obj=nl.LMA2.0R, trace.progress=trace.progress,
                              parameter.names=input.names, iter.length=iter.length,
                              function.name="Morris")}
                           )

# iterate over evalulation criteria

par(mfrow=c(7,2),  mar=c(4,4,2,1))
titles = c("unemployment density", "tightness fluctuation", "ltu density", "wages", "ltu wages", "duration", "ltu duration")
for (i in (1:7))
{
  tell(mo, sim.results.morris[i,])
  mu <- apply(mo$ee, 2, mean)
  mu.star <- apply(mo$ee, 2, function(x) mean(abs(x)))
  sigma <- apply(mo$ee, 2, sd)
  
  # mu* against mu
  if (plot.on.screen == TRUE) {
    par(ask=TRUE)
  }
  plot(mu.star, mu, pch=1:6, col=1:6, axes=F, ann=F, cex=4, lwd=5)
  box()
  axis(1, cex.axis=2)
  axis(2, cex.axis=2)
  title(main=titles[i], font.main=3, cex.main=2)
  title(xlab='Morris mean value (mu*)', cex.lab=2, line=2.1)
  title(ylab='mean value (mu)', cex.lab=2, line=2.3)
  
  # mu* against sigma
  if (plot.on.screen == TRUE) {
    par(ask=TRUE)
  }
  plot(mu.star, sigma, pch=1:6, col=1:6, axes=F, ann=F, cex=4, lwd=5)
  box()
  axis(1, cex.axis=2)
  axis(2, cex.axis=2)
  title(main=titles[i], font.main=3, cex.main=2)
  title(xlab='Morris mean value (mu*)', cex.lab=2, line=2.1)
  title(ylab='sigma', cex.lab=2, line=2.3)
}